Registration of Tungsten fibres on XCT images

This demo aims to demonstrate the use of gVirtualXRay and mathematical optimisation to register polygon meshes onto X-ray microtomography (micro-CT) scans of a tungsten fibre. Below is an example of CT slice.

The fibre.

Our simulations include beam-hardening due to polychromatism and they take into account the response of the detector.

We use SimpleGVXR's Python wrapper and Python packages commonly used in tomography reconstruction (Tomopy), image processing (scikit-image and SimpleITK), computer vision (OpenCV), and non-linear numerical optimization (CMA-ES, Covariance Matrix Adaptation Evolution Strategy).

Registration steps

  1. Initialisation
  2. Registration of a cube
  3. Optimisation of the cores and fibres radii
  4. Recentre each core/fibre
  5. Optimisation the radii after recentring
  6. Optimisation of the beam spectrum
  7. Optimisation of the Poisson noise
  8. Optimisation of the phase contrast and the radii
  9. Optimisation of the phase contrast and the LSF

Import packages

We need to import a few libraries (called packages in Python). We use:

Global variables

We need some global variables.

Load the image data

Load and display the reference projections from a raw binary file, i.e. the target of the registration.

In the literature, a projection is often modelled as follows:

$$P = \ln\left(\frac{I_0}{I}\right) = -\ln\left(\frac{I}{I_0}\right) = \sum_n \mu(n) \Delta_x$$

reference_normalised_projections loaded from the binary file corresponds to $\frac{I}{I_0}$. The flat-field correction has already been performed. It is now necessary to linearise the transmission tomography data using:

$$-\ln(normalised\_projections)$$

This new image corresponds to the Radon transform, known as sinogram, of the scanned object in these experimental conditions. Once this is done, we divide the pixels of the sinogram by $\Delta_x$, which is egal to the spacing between two successive pixels along the horizontal axis.

We define a new function to compute the sinogram from flat-field correction and calls it straightaway.

CT reconstruction

Now we got a sinogram, we can reconstruct the CT slice. As we used a synchrotron, we can assume we have a parallel source. It means we can use a FBP rather than the FDK algorithm.

Normalise the image data

Zero-mean unit-variance normalisation is applied to use the reference images in objective functions and perform the registration. Note that it is called standardisation (Z-score Normalisation) in machine learning. It is computed as follows:

$$I' = \frac{I - \bar{I}}{\sigma}$$

Where $I'$ is the image after the original image $I$ has been normalised, $\bar{I}$ is the average pixel value of $I$, and $\sigma$ is its standard deviation.

Set the X-ray simulation environment

First we create an OpenGL context, here using EGL, i.e. no window.

We set the parameters of the X-ray detector (flat pannel), e.g. number of pixels, pixel, spacing, position and orientation:

3D scene to be simulated using gVirtualXray

And the source parameters (beam shape, source position)

The beam spectrum. Here we have a polychromatic beam, with 97% of the photons at 33 keV, 2% at 66 keV and 1% at 99 keV.

The material properties (chemical composition and density)

The LSF

Find circles to identify the centre of fibres

We can use the Hoguh transform to detect where circles are in the image. However, the input image in OpenCV's function must be in UINT8. We blur it using a bilateral filter (an edge-preserving smoothing filter).

Convert the image to UINT8

We first create a function to convert images in floating point numbers into UINT8.

We blur the CT scan using a bilateral filter. It preserves edges.

Apply the Hough transform

Overlay the detected circles on the top of the image

Unlike the previous example, did did not work that well. Here 13 fibres were missed. Many centres are also misplaced. We will use another technique to register the fibres, the popular Otsu's method. It creates a histogram and uses a heuristic to determine a threshold value.

Clean up

Mark each potential tungsten corewith unique label

Object Analysis

Once we have the segmented objects we look at their shapes and the intensity distributions inside the objects. For each labelled tungsten core, we extract the centroid. Note that sizes and positions are given in millimetres.

We now have a list of the centres of all the fibres that can be used as a parameter of the function below to create the cylinders corresponding to the cores and the fibres. For each core, a cylinder is creatd and translated:

gvxr.emptyMesh("core_"  + str(i));
        gvxr.makeCylinder("core_"  + str(i), number_of_sectors, 815.0,  core_radius, "micrometer");
        gvxr.translateNode("core_"  + str(i), y, 0.0, x, "micrometer");

For each fibre, another cylinder is created and translated:

gvxr.emptyMesh("fibre_"  + str(i));
        gvxr.makeCylinder("fibre_"  + str(i), number_of_sectors, 815.0,  fibre_radius, "micrometer");
        gvxr.translateNode("fibre_"  + str(i), y, 0.0, x, "micrometer");

The fibre's cylinder is hollowed to make space for its core:

gvxr.subtractMesh("fibre_" + str(i), "core_" + str(i));

Registration of a cube

We define a function to create the polygon mesh of the Ti90Al6V4 matrix.

Simulate the CT acquisition

  1. Set the fibre and cores geometries and material properties (Step 39)
  2. Set the matrix geometry and material properties (Step 40)
  3. Simulate the raw projections for each angle:
    • Without phase contrast (Line 5 of Step 45), or
    • With phase contrast (Lines 14-55 of Step 45)
  4. Apply the LSF (Lines 57-60 of Step 45)
  5. Apply the flat-field correction (Step 62)
  6. Add Poison noise (Step~\ref{??})
  7. Apply the minus log normalisation to compute the sinogram (Step 63)

Compute the raw projections and save the data. For this purpose, we define a new function.

Flat-filed correction

Because the data suffers from a fixed-pattern noise in X-ray imaging in actual experiments, it is necessary to perform the flat-field correction of the raw projections using:

$$corrected\_projections = \frac{raw\_projections\_in\_keV − dark\_field\_image}{flat\_field\_image − dark\_field\_image}$$

Note that in our example, $raw\_projections\_in\_keV$, $flat\_field\_image$ and $dark\_field\_image$ are in keV whereas $corrected\_projections$ does not have any unit:

$$0 \leq raw\_projections\_in\_keV \leq \sum_E N_0(E) \times E\\0 \leq corrected\_projections \leq 1$$

We define a new function to compute the flat-field correction.

The function below is used to simulate a sinogram acquisition. Phase contrast in the projections can be taken into account or not.

The function below is used quantify the differences between two images. It is used in the objective function.

The function below is the objective function used to register the matrix.

Apply the result of the registration

Simulate the correspond CT acquisition

Display the result of the registration as an animation

Animation of the registration (GIF file)

Adding the fibres

The radius of a tungsten core is 30 / 2 um. The pixel spacing is 1.9 um. The radius in number of pixels is $15/1.9 \approx 7.89$. The area of a core is $(15/1.9)^2 \pi \approx 196$ pixels.

Optimisation of the cores and fibres radii

The function below is the objective function used to optimise the radii of the cores and fibres.

Animation of the registration (GIF file)

Apply the result of the registration

The 3D view of the registration looks like:

3D view

Recentre each core/fibre

Each fibre is extracted from both the reference CT slice and simulated CT slice. The displacement between the corresponding fibres is computed to maximise the ZNCC between the two. The centre of the fibre is then adjusted accordingly.

Applying the result of recentring

Optimisation the radii after recentring

After recentring the centres, another run of optimisation is executed to refine the radii of the fibres and cores.

Animation of the registration (GIF file)

Apply the result of the registration

Optimisation of the Poisson noise

Apply the result of the optimisation

Optimisation of the beam spectrum

Animation of the registration (GIF file)

Optimisation of the phase contrast and the radii

Animation of the registration (GIF file)

Apply the result of the registration

Optimisation of the phase contrast and the LSF

Animation of the registration (GIF file)

Apply the result of the registration

Extract the fibre in the centre of the CT slices